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ABSTRACT 


Title: PRESSURE FIELD DUE TO UNSTEADY LOADING OF A MARINE PROPELLER 


Author: J.ReP. D'Almeida 


Submitted to the Department of Naval Architecture and Marine 
Engineering on May ik, 1971, in partial fulfillment of the require- 
ments for the degrees of Naval Engineer and Master of Science in 
Naval Architecture and Marine Engineering. 


A vortex model is used to develop a scheme to compute the 
pressure field around a marine propeller operating in a wake which 
is nonuniform in both the radial and circumferential direcel Oils. 
This model takes into account the effects of thickness form and 
outline of the blades and the form of the unsteady loading distri- 
bution, both radially and over the chord. The model can be used 
for lightly as well as "moderately" loaded propellers. Following 
the theoretical development a computer program is written in 
FORTRAN IV, to be used with the IBM System/360-65 digital computer. 


The workability of the program 1S tested with sample calculations. 
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Chapter I 


INTRODUCTION 


The propeller induced vibratory forces acting on a ship have 
been for a long time an important subject of experimental and theo- 
retical research. Although some important conclusions have been 
reached based on reliable experimental results ie it is always 
hard, if not impossible, to support such conclusions with theoretical 
calculations. One reason for this defficiency is that for a long time 
no analytical method was available for evaluating the propeller loading 
as a three-dimensional problem. 

The initial numerical results obtained for the propeller-induced 
pressure field were based on a simple lifting-line representation [10] 
which was found to predict pressures far below those obtained from 
measurement. This theory was later refined by including the effect of 
the blade thickness and this brought the numerical results closer to 
the experimental data [9, 10, 11]. 

The effects on the pressure field due to nonuniform inflow were 
studied by Tsakonas et al. [8]. They treated the problem as a Dirichlet 
problem of potential theory and related the pressure field near the 
propeller to the trust- and torque-producing forces. These forces 
were computed using the gust theory for airfoil, and assuming that 
the flow around each blade was two-dimensional. Two correction factors 
were added to compensate for cascade and finite aspect ratio effects. 
Even though the calculations made using this method cannot be con- 


sidered very accurate, they can give an idea of the relative importance 
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of the unsteady effects arising from the nonuniformity of oe ship's 
wake. This can be seen in Fig. 1, taken from [8]. It is apparent 
that very close to the plane of the propeller the increase in pressure 
due to nonuniformity is at most 20 percent of the total pressure but 
beyond x/D = 0.4 the percentage change becomes much greater. The 
reason is that at such distances all the pressures are small, but the 
slow rate of decay of the components arising from nonuniformity makes 
them more important. It must be noted, however, that the pressure due 
to thickness is not included in Fig. 1! so that the relative effect 
close to the propeller arising from nonuniform inflow is of the order 
of 10 percent for the particular wake used in these calculations. 

A truly three-dimensional theory was recently used by Kerwin 
[18] to compute the field point velocities induced by a marine pro- 
peller operating in a steady flow. He makes use of the vortex theory 
to represent the blade loading and derives a scheme to compute the per- 
turbation velocities on specified field points located on the blades. 

In the present work we will extend the program developed by 
Professor Kerwin to include field points located at any point in the 
fluid around the propeller and to take into account the unsteady effects 
due to a nonuniform inflow such as is found in a ship's wake. In 
addition, a scheme will be developed to compute the velocity potential 
at the field points. The pressure field is then computed by sub- 
stituting the perturbation velocities and time variation of the 


velocity potential into Bernoulli's equation. 
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The usual assumptions of an incompressible, frictionless and 
irrotational flow are held. It is also assumed that all the disturb- 
ances are small, which requires that the blades be thin and propeller 
loading small. The effects of thickness form and outline of the 
blades and the form of the loading distribution, both radially and 
Over the chord are taken into account by this theory. 

When this work was already under way, the Davidson Laboratory 
published a report by Jacobs et al. [19] which also treates the un- 
steady, propeller induced, pressure field weine a lifting surface 
theory. However, their treatment is based on the solution of a Dirichlet 
problem for the pressure field. It seemed of interest to compare the 
results obtained with the two different approaches but unfortunately 
were were not able to obtain from D.L. the loading distribution for 


the propellers used in their report. 





Chapter II 


FORMULATION OF THE PROBLEM 


It is desired to find the pressure at any point in the fluid 
around a ship's propeller. Assuming that the flow is incompressible, 
frictionless and irrotational, the problem becomes one of potential 
theory which in principle can be solved by finding the velocity 
potential and velocities induced by the propeller and substituting 


into the Bernoulli's equation to solve for the pressure. 


lL. Géemetry sof the Propeller 


The following geometrical considerations are adopted from [1]. 

Two orthogonal coordinate systems -- a Cartesian (x,y,z) and a 
cylindrical (x,r,9@) -- are associated with the propeller as represented 
in Fig. 2. Both systems are fixed with the propeller and the first 
blade is made to coincide with the y axis. All length dimensions are 
nondimensionalized with respect to the radius of the propeller, R. 


The equations relating the two coordinate systems are: 


a 
2 


(y” <5 2°) 


r cos@ 1s 


K< 
I 


(1) 


r sind 6 tan (z/y) 


N 
Il 


In order to relate corresponding points on each of the K blades, 
we define oy as the @ coordinate of the point at the tip of the k'th 


blade. For symmetrically arranged blades it is 


oe an (2) 
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FIG. 3 = PROJECTED BLADE OUTLINE 
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This is illustrated in Fig. 3 which also shows the projected blade 
outline, defined by the angular coordinates 6, (x) and 6,,(r) of the 
leading and trailing edges, respectively. 

The geometry of the flow past a blade is schematically re- 
presented in Fig. 4. It is seen that the oncoming flow V,@ forms 


an angle 8 with the yz plane, and we define 


Nor) = tanb(r) (3) 
as the "'advance coefficient''. The streamlines are helices with 
pitch 

PCa) = 27 AG) (4) 


It is useful to define an orthogonal curvilinear coordinate 


system (s,r,n) associated with the helical surfaces swept out by 


the undisturbed flow past the radial lines {x = 0, @ = 5}. This 
system, shown in Fig. 2, has unit vectors: 

i 7 sinB - sin(@ +8, )cosB +k cos(@ +8, )cosPB 

aL = cos(@ +5,) +k sin(@ +8,) (5) 


a cosB + sin(6 +6,)sinp =e cos(@ +5, )sinB 


It is assumed that the blades lie on those helical surfaces. 
This assumption is in accordance with the linear theory for lightly 
loaded propellers, however, it does not impose any limitation to the 
possibility of including nonlinear refinements by suitably modifying 


the advance angle [1]. 

















FIG. 4 - APROACH FLOW VELOCITY DIAGRAM 





~10- 


2. Expression for the Pressure 


The Bernoulli's equation for unsteady flow is 


p=, 5 +40,(, - V) 
where 
p = hydrodynamic pressure 
ae fluid density 
® = velocity potential 
Va = speed of free stream or speed of advance 
V = resultant velocity at a point 
t = time variable. 


Changing to a frame of reference moving with the propeller 
blades and after linearization the above equation becomes [11] 
op 
Des D- e + p uv, + p Or u. (6) 
where 


axial and tangential disturbance velocities 


an 
I 


re) 
ll 


angular velocity of the propeller. 


3. Induced Velocities 

In order to find the velocities induced by the propeller, we 
make use of the vortex representation described, for instance, in 
[1] and [2]. The following is an adaptation to the unsteady case 


of the relations derived in [1]. 
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The inflow velocity and verticity strength are expressed as 
complex Fourier series of the angular position of the rotating co- 


ordinate system with respect to a set of axes fixed in the ship, 


V Co pt) ~ TGs am 
r(o,0,t) = 7(p,0)-e 
where 
Ue = Seaeeaine flow velocity past the blade section 
Y = vorticity strength, in general 
N = harmonic number 
Q = angular velocity of the propeller 


(i) Velocity Induced by Source Distribution 
The velocity induced at a point (x,r,@) by sources distributed 


over K blades is 


(s) Kj iN(2t-8x) Y ee ' 
u, Gor .e,t) =2 ia | \ V 69) 
k=l Ty Oo, 


Pie ena ec en (a) 1 
xe iN — _ {x (po) 5 Nee) | dp do (7) 


iN(Qt-5,) a) 
(s) K a i\ T z 
u Gaeoet) =22 — wae \ ve) 
a Ty 91, 


i 
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ys. 6 Sint, -6)-[p + (0) 


Sea Cm OU FDU 


iN, 





Sie 


where 


_ 
2 


D 


[ (x -A(p) moe wes a 2rp cos(@+ oe 6) ] 


blade thickness 


5 


The superscript s denotes a source velocity. 


(ii) Velocity Induced by Vortex Distribution 


We separate the components of the induced velocities due to 
the bound vorticity in the blade (superscript b), radial vorticity 
in the wave (sup. r), streamwise vorticity in the blade (sup. f) 


and streamwise vorticity in the wake (sup. w). 


On : 
\ V1, (9) i 


fe “1. 


uO) (x, 2,0 ,t) =) = 


“TT 
k=1 


a — 
tip + A (po) ] sin@+68, - @) 


2 ae dp dp (9) 
D 
(b) K iN (28 ,) z on 
a CoE Gee) =) Sa { \ Y(0 59) - 
Ser Th OL 


a 
2 


(x- Mpolo + *(o) }* cos(+8, - 6) 


ee ee A EC) 
D 
, 1 © 
K  iIN(Qt+OT-d,) ¢ : 
my 8,0) =), = ‘i \ \5 (o) 20 
a el if r 
Yn ST 


54 
2 


ee No) | sin(o+,- 8) 
: aa dp do (11) 
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(x- Mop) [p+ 7) ]* cos(+8,,- 8) 


3 dp do (12) 
D 
IN (2Qt-8 
0 .& ( ye Op 
(x,r,0,t) =2) =) an Y (0,9) 
ead 
Th OL 


Bes - pr cos(pt Sis Q) | 


a dp do (13) 
K gREESS a 
ey 23038) 2 — St aa 
Th OL 


{ Mp) [r- p cos(pt+ oe Q) | 
° re of 


3 
D 


o[x- A(p)p] sin(ot+ 8, - 6) 
+ ———5*—} ap ap (14) 
D 
: IN(Qt+0,.-5, ) 
a) (x,x,6,t) = —— ie (0) °p ae 
k=] r. 6 
aia 


Z 
[p - pr cos(@+8,- 6)] 


3 dp dp (15) 
D 
IN(Ot+9,,-5) a 
a) (x,2,8, £) = 2) os ie i \ io) a we 
k=l rico 


Mp)[r- p cos(pt ay @) of{x- Mo)p}sin@+s, - @) 
D D 
(16) 





ihe 


where 


¥, (p.9) = complex amplitude of the bound vorticity 
¥,(p) = complex amplitude of the radical vorticity 
in the wake 
¥ (050) = complex amplitude of the streamwise vorticity 
in the blade 
¥_(o) = complex amplitude of the streamwise vorticity 


in the wake. 


3. Velocity Potential 

As before we separate the effects due to source and vortex 
distributions. 

(i) Velocity Potential due to Source Distribution 

The velocity potential at a point due to a source of unit 


strength located at a distance, d, is given by [16] 


Multiplying by the source strength, integrating over the first 


blade and summing over the K blades gives: 


le, 


~ in(%t-8, ) 1 5 
i Geen), 
i \ V.,60) *e aS 


(s) -e 
@"’ (x,r,6,t) =2, —~———_ 
a War 


2 2 3 
oes 2 (o) |] Py Glo (17) 
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(ii) Velocity Potential due to Vortex Distribution 

Instead of trying to find an exact expression for the velocity 
potential, which would be rather cumbersome, we make use of the fact 
that the unsteady effects become important only at a distance from 
the propeller and that a reentrant vortex is, in a sense, equivalent 
to a uniform distribution of double sources over any surface bounded 
by it. It is shown in Appendix A that if we substitute a reentrant 
vortex by an equivalent concentrated doublet, the error in the value 
of the velocity potential decreases rapidly with distance. 

In this way we can compute the velocity potential by summing 


a set of discrete values given by expressions of the form (see Ap- 


pendix A) 
one wa: e4 
T DH 
where 
S = strength of the concentrated doublet 
Y = angle between the doublet axis and the line from 


the doublet to the field point. 
The value of S will be found later and cos ¥ is easily computed 


from 


ao | 
or 





cos Y = 


SI 
on 


where 


vector direction of the double axis 


5 
I 


= 
I 


[x- A(p)p]it [xr cosO- 0 cos(p +8.) ]j + 


+ [r sind- » sin@+ 5.) ]k 





SG. 


Using (5) to obtain n, it follows that 


[x - \Co)p JcosB(p) +r sinB(p)cos(p+8, - @) 
|, (19) 


4. The Loading Modes 

For the purposes of this work it is assumed that the distri- 
bution of the loading on the propeller blades has already been de- 
termined by any of the known methods (examples [2], [3], [7], [12] 
and [15]) which solve the unsteady lifting surface problem of a 
propeller operating in a wake. The usual procedure in these methods 
is to assume that the loading may be expressed in separable form as 
the product of a function of the spanwise coordinates and a function 
of the chordwise coordinates. In selecting these functions, one 
makes use of the known results of the theory for incompressible flow 
in two dimensions and of the lifting line theory for the higher 


aspect ratio wings. This leads to functions of the type 


[9 @) 
ae 2; FCs) x G(r) 
n=O 
where 
Ap = loading per unit area 
FCs) = chordwise function 
Gr) = spanwise function. 


Two dimensional airfoil theory suggests that we choose Fe on 


the form, [5] 





F C8) = E > j + sin[n(cos’—)] 





ee 


where 


2(s - Sie) _ 4 
L 


£ 


chord length 


Introducing a new coordinate 6 defined by 


S 7 Sie £ A 
$s =z OF C088 0 


LA 
D> 

LA 
=} 


Ee becomes 


“A 


F_(6) = @etan S _ cos(nir) sin(nd) 


The spanwise function G(r) is either determined directly with 
the aid of a specially selected interpolation function [5] or ex- 


pressed in the form 


ae A ae (r) 
m 
and then determined by the solution of the integral equation resulting 
from the formulation of the lifting surface problem. [In order to ob- 
tain a numerical solution Gf) is usually approximated by a stepwise 
distribution. 


The loading then will be of the form 


NID») 


+D L(x) sin(nd) (20) 
n=] 


Ap(r ,0) = L@) cot 
where L, and La are stepwise functions of the radius. This is called 
the Birnbaum distribution. 


When the loading is expressed in terms of the bound vorticity, 


it is necessary to add an extra term to the Birnbaum modes in order 
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to satisfy the Kutta condition at the trailing edge [13]. In this 


case it becomes 


Nol@ > 


+D A (rx) sin(nd) + £(r,6) Cay 


¥,, (758) = A, (x) cot 
n=1 


where the Kutta condition requires that 


f(r ,7) = 7, (1,7) = “ik I (22) 
where 

, = N& LO) 

n 2V (x) 

I’ = chordwise integration of ¥ 


b 


According to Brown [3] a convenient choice for f(r ,@) is 


1 - cose 
2 


£(r ,0) = -ik Ix 

After substituting into (21) and integrating over the chord, it fol- 
lows that 

ik 

n 


T+ ik or 
19) 


a { - cosé 
x (A, +2 A,) ea a 


A 
f(r,6) = -T x 
It must be noted that there is no correspondence (except for 
the zeroth harmonic) between the coefficients of (20) and (21). 
However, when the loading is given in terms of pressure distribution 


one can always find the corresponding vorticity distribution from 


the solution of the following Volterra equation of the second kind 


[14] 
5 ° 
oi aaa 2) Y,(r,s8,t)ds + V(r) ¥,(r,8,t) (24) 
t 


aA 
in terms of @ the above expression becomes 





mii 


ray 


er Ap(r ,@) 
v7 59) . P, Vi) 


0 
: ik (x) Vietd V,(©,8) de (25) 
0 
5. The Vorticity Relations 
The pulsating bound vorticity is expressed in the form 
Y, 60.9,k, t) ia Y, 60) e (26) 
All the other components of the vorticity can be expressed in 
terms of ee as required to preserve continuity of vorticity (see [1] 


and [2] for detailed derivation) 


iN(Q2t-5,) 
¥,(0,9,k,t) aes ¥,C0,9) S (27a) 
7 -(0,0) = \" iG tp? + Moe) 7? 
¥_{0,9) = - 5 oP ¥,(0:P)*[o + A Cp) }* }dp (27b) 
L 
iN[ Mt -@-9,,) -5, J 
Y,,(0,9,k, t) = Yo) e (28a) 


¥_,(o) = - oe (28b) 


iN[2t -(p-6,,) -8, 1 


¥(©,9,k,t) = r,(p) (29a) 
¥,(p) Bee (o)e (29b) 


Pa 2 4 
on cp) | 
where 


k. T_ 7 9 a 
To) = \ 40,9) *fo°+ M(@)]? ap (30) 
OF 





~-20— 


Notice that an and ie (streamwise and radial vorticity in the wake) 
are expressed as the product of a complex amplitude at the trailing 
edge and a sinusoidal time and distance variation along the helical 
surfaces in the streamwise direction. This form had already been 


implied in eqs. (6), (7), (10) and (11). 
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Ciapeer Til 
NUMERICAL SOLUTION 


The following is an outline of the computational model used 
to solve the equations derived in the previous chapter. 

The blades and wake are divided into chordwise panels and 
these panels are Abe aed by radial grid lines as shown in Fig. 6. 


The scheme for dividing the span is shown below. 


ad/fh da/2 d 


FIG. 5 - SPAN SUBDIVISION FOR NUMERICAL SOLUTION 


The interval is first divided into MN equal spaces. The intervals 
at each end are then subdivided into half and quarter spaces as 
shown. The result is MT = MN+4 chordwise panels. The chord length, 
Pucca weeahe. | bound ciretwlation strength and source strength are as- 
sumed to be independent of radius within each of the panels, with 
values matching those of the continuous functions at the mid-radius. 
The continuous trailing vortex sheet then becomes a set of MI+ 1 


concentrated trailing vortex lines at each of the panel boundaries. 
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FIG. 6 = BLADE AND WAKE LATTICE FOR NUMERICAL SOLUTION 
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The bound circulation in the blade and the radial circulation 
in the wake are assumed concentrated ateach of the radial grid 
laine s. 

The above assumptions allow us to write the integrals (7) to 


(12) as summations, as follows: 


: MT I ; ; P 
-_ eee) = iNet > acip) --a + [x - A(m) -p(i) ] 
1= 


m=1 1 
K -iNS mt+l 
k d 
De ( 5 (31) 
=| 'e D 
m 


; Mi i ore ee ; 
uf) (x,r,8,¢) _ pi NOt sy > 27 iNp (i) SBtn 4) 








no) 
K -iN5, a 
2) e wo eniaiG). + Syme?) le) = (32) 
=F - D 
jaal 
MT I 
0) (x,v,0,t) = gee 2 a ae 2) 
m=l1 i=l 
K -iN5, Pmt] 7 
2 e sin(p(i)+58, - 6) \ =a (33) 
Pe) k D> 
Pn 
MT I 
b INQ 
we (x,r,0,0) = -e NEF SBOLY oy  AGm) pCi) 
m=l1 i=1 
2 Pmt1 
4 d 
2 6 -cos(p(i)+3, - @) 4 —- (34) 
k=1 k p> 
Pm 
° MT oc iN (6 ay) (1), 
Me Cenn0.,t) — _ Qi Nt “re >> yp e tT ° GR(m) 
= m=L i=] 2 


K ~ iN5 |. ce do 
Qu e -sin(@p(i) +8, - 0) - 3 (35) 
k=1 on D 
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: MT co iN(O -p(i)) 
INQt iE GR 
vw (x,r,0,t) = -e NTO De , ROD Gx - Mm) «9 (4)) 
m=] i=I1 
Dai a Goe(OCi)ytG.=0) - { sale (36) 
= k 3 
k=1 O D 
m 
where 
SB(m,i) = amplitude of the equivalent source strength, 
concentrated at the i'th grid line in the m'th panel. 
GB(m,i) = amplitude of the equivalent bound vorticity, con- 
centrated at the i'th grid line in the m'th panel. 
GR(m) = amplitude of the radial vorticity in the wake, con- 


centrated at the i'th grid line in the m'th panel. 
In the above equations the source and vortex strengths are non- 
dimensionalized with respect to 2TRV . and the perturbation velocities 
are non-dimensionalized with respect to ee 
The integrals that appear in the equations above are evaluated 
between the limits of integration corresponding to the boundaries 


of each panel and can be determined analytically as follows 


Pmt] O 
: dem rp cos( (i) + 8,- 6) + (x- A(m) eC) + o ie 
o, DY ([x- Mm)o(i)] + r7sin?@(i) +8, - @))-D 
m 
Pit] : le) 
{ ds po - x cos(p(i) se One @) m+] 
cas mo? a ave 
p, dD [(x - A(m)p(i))~ + rv sin” @(i) +5, -6)]-D 
m 
The contribution from the trailing vortices -- eqs. (13)-(16)-- 


raises the problem of streamwise integrals that cannot be evaluated 
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analytically. It was found by Kerwin [18] that for the zeroth 
harmonic a simple trapezoidal integration is as efficient as more 
sophisticated quadrature formulas, except for nearby elements. In 
order to investigate the application of this argument to other har- 
momics, we carried out the integration for a single trailing vortex 
line located at 0./75R, using a Gaussian quadrature and a trapezoidal 
rule. The computation of the perturbation velocities was carried 

for different points around the propeller using a third harmonic of 
the loading. We found that the results obtained with the two for- 
mulas differ in all cases by less than 2 percent. However, due to 
the pulsating character of the wake vorticity it is necessary to keep 
the intervals of integration small compared with the wavelength of 
the perturbation. This results in a very large number of computations 
for each fortex line (about 1100 for a field point in the plane of 
the propeller). A better alternative was found by substituting the 
infinite integrals by an infinite series of finite integrals and then 
applying the Euler transformation to speed up the convergency of the 
series. Basicly the method works as follows: 


Given an infinite integral of the type 


1¢.@) 


I = if Soaks *-f(x) dx 
0 


t—Momtranstonmed into the infinite Series 


(k+1) 1/N 
as Bee f(x) ds = 
k=0 “kK a/N 
Bs 1/N 
= JE, er ore + k =) dt 
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Now, applying the Euler Transformation [21], the above series 


becomes 
90 (aon a, k lene 
I — ST ; Aa, —- (=1) ) a 
1/N 
a \ NG sa gerea T/N) dt 


Sample computations done with this method showed good convergence 
with a fairly reduced number of terms. Also, it was found that a 
5-point Gaussian quadrature was adequate for the computation of the 
tenpms of the series. 

When the above transformations are applied to eqs. (15) and 


(16), we obtain 
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- Gaye a 
: ee (38) 
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-iN(E+0,,) 
ay = e FWI(E+ 6,4 k /N) dé 
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where 
TL(m) = amplitude of the m'th trailing vortex, on the wake 


ofa Efe er r cos(gt+ Sie 6) | 


D> 


It 


FWA (9) 


FWI@) = {-p [A(m)cos(p+ 6, - 6) - (x- A(m) -—) sin@@+5, - @)] + 


EEN De 


It is clear that eqs. (37) and (38) cannot be applied to the zeroth 
harmonic of the loading because the upper limit of the integrals 
would go to infinity. Instead, we will use the trapezoidal rule 
which was found appropriate for this case, as mentioned before. 


Using a trapezoidal rule, eqs. (15) and (16) become 


MT+1 ©  iN(O,-p(i)) : 
oY, cre = eiNOt 3 a T ; Re 1) 
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~iN5 ~r cos((i)+ By @) 
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eat D 
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ae Cle z oi Nt | > Sc we 
m=1 i=I 
K ~iNS, 
piace (-e [Atm cos(p(i) +8, - @) - 


(40) 
The contribution from the trailing vortices on the blade -- eqs. (13) 
and (14) -- are computed using a trapezoidal integration, for all the 


harmonics, 
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cE e 
m={ 1i=1 


K ~1N5,. 
> ere {-p [A(m) cos(p (i) + 8, - Q) - 
k= = 


—» 


(x - (mp (4) *sin(@@(4)+8,- O)] +r Mm))/D™ (42) 
where 
GT(m,i) = Amplitude of the m'th trailing vortex located 
on the blade, between the i'th and the i'tht 
erid lines. 


The triple summations corresponding to the velocit otential 
P P & y Pp 





due to the source and vortex distributions -- eqs. (17) and (18)-- are 
, MT I : : ; 
Dea eo) _ Qt Nat iy hee inp (i). ee : 
m={ i=1 
K  -iNS » Prt] 
k dp 
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D 
(45) 
where 

a” = velocity potential due to the vorticity in the blade. 

o = velocity potential due to the vorticity in the wake. 
GP(m,i) = amplitude of the concentrate vortex equivalent to 

GB(m,i). 

GL(m,i) = amplitude of the concentrate vortex equivalent to 


GR(m,i). 
The integral appearing in eq. (43) is evaluated analytically 


mee Prckl 
= oh [22+ p- x4 cos(p(i) +8, - 6))| 
5 p 


m m 
The value of the bound circulation is obtained either directly 
if the loading is expressed as in eq. (21) or after solving a Volterra 
equation (25) if it is expressed as in eq. (20). 
The solution of the Volterra equation can be obtained using the 


Liouville-Neumann method of successive substitutions [20]. This gives 


pelle temo eyo 
Tees eR. = 5p Vi) [A pr ,8) 7 GE ©) ~- 


© (ak (x) 87 —  cik(r) ¢? 
pe -i) - ae cos@ e © J Ap(r ,t)sintdt] 
k=} 0 

(46) 
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When the pressure loading is given in the form of Birnbaum modes 


(20), the above integral is evaluated using the relations 


6 
i cotan ie sin t dt = 6 - aoe) 
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Chapter LV 


RESULTS 

The numerical solution developed in the preceding chapter was 
Puc into computer form; the details of the program are presented in 
Appendix B. In order to test the workability of this program we 
used it to compute the pressure at an upstream field point using the 
NSRDC Propeller No. 4118 [22]. The design characteristics and oper- 
ating conditions of this propeller are represented in Figs. 7 and 8. 
We considered both the open water condition and the nonuniform wake 
obtained with a screen which produces a circumferentially varying 
sinusoidal wake composed only of multiplies of blade harmonics. 
Unfortunately, we did not have information about the loading distri- 
bution for those conditions and therefore we used an arbitrary number 
of modes of the Birnbaum distribution with coefficients equal to 
fart. In order to test the solution of the Volterra equation (46) 
the input loading was expressed in terms of pressure distribution 
rather than vorticity distribution. | 

We made several sample calculations in order to evaluate the 

\ 

effect of changing the grid spacings and to determine how far down- 
stream we needed to carry the streamwise integrations in the wake 
(Table 1 shows a sample of the computer output for one of the cal- 
culations). We found that for the field point considered (x = -.8, 


r = 1.2, 8 = 0), we needed about 20 chordwise panels and had to keep 


the angular grid spacing down to 2 degrees. The contribution due to 
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FIG. 7 ~- PROPELLER 4118 DESIGN CHARACTERISTICS 


1/D 
0.3200 
0.3635 
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0.4392 
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op on wr ep aw os 


T/1 
02026 
O707'7'6 
0.0590 
0.0451 
0.0348 
0.0271 
0.0210 
0.0166 


P/D 
1.0860 
1.0845 
1.0823 
1.0796 
1.0780 
1.0766 
1.0750 
1.0735 
1.0700 
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FIG. 8 = RADIAL DISTRIBUTION OF THE AMPLITUDES AND 


PHASES OF THE THIRD HARMONIC OF THE 


CIRCUNFERENTIAL WAKE VARIATION 
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“ewer == PRESSURE FIELD AROUND A MARINE PROPELLER **#os sexe 


PROPELLER CHARACTERISTICS 


ee, NSRDC, PROPELLER NO.4118 


DIAMETER 12.0000 INCHES 
HUB DIAMETER 2.4000 INCHES 
NO. OF BLADES 3 


RPM 1080.0000 


BLADE OUTLINE (UNITS-FRACTION OF RADIUS) 


RADII LEAD.EDGE TRAIL.EDGE MAX.THICK. 
64250 -0.3417 0.3417 0.0605 
0.300 -0.3635 0.3635 0.0564 
0.400 -0.4048 0.4048 0.0478 
0.500 -0.4392 0.4392 0.0396 
0.600 -0,4610 0.4610 0.0321 
0.700 -O.4615 00,4615 0.0250 
0.800 -0.4347 0.4347 0.0183 
0.900 -0.3613 0.3613 0.0120 
0.950 -0.2500 0.2500 0.0065 
1.000 0.0 0.0 0.0 


GRID SPACING 


NO. OF CHORDWISE PANELS 20 
ANGULAR CHORDWISE SPACING 0.0349066 RADIANS 


TABLE 1 - SAMPLE COMPUTER OUTPUT 


a es 


INFLOW CHARACTERISTICS 


SYMBOLS AND UNITS 


VS SHIP OFPEED, SeNOTS 
RADII FRACTION OF RADIUS 
AVXO ZEROTH HARMONIC OF INFLOW / VS 


TANBETA TANGENT OF ADVANCE ANGLE 

NHARM HARMONIC NUMBER 

XVAN AMPLITUDE OF NHARM OF INFLOW / VS 
XPHN PHASE OF NHARM OF INFLOW, RADIANS 


VS= 8.8780 NHARM= 3 


RADII XVXO -- TANBETA XVXN XPHN 
0.250 1.0000 1.0606 0.0800 ~-0.2269 
0.300 1.0000 0.8838 0.1000 ~-0,1920 
0.400 1.0000 0.6629 0.1500 -0.1047 
'eesa0” 1.0000 “0.5303  omee00 0.0 
0.600 1.0000 0.4419 0.2500 0.0873 
0.700 1.0000 0.3788 0.2900 0.0 
0.800 1.0000 0.3324 0.3200 -0.0324 
0.900 1.0000 0.2946 0.3400 0.0 
0.950 1.0000 0.2791 0.3450 0.0436 


TABLE 1 (cort.) 





RADIT 


0.250 
0.250 
0.250 
0.300 
0.300 
0.300 
0.400 
0.400 
0.400 
0.500 
0.500 
0.500 
0.600 
0.600 
0.600 
0.700 
0.700 
0.700 
0.800 
0.800 
0.800 
0.900 
0.900 
0.900 
0.950 
0.950 
0.950 
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LOADING DISTRIBUTION 


(VORTICITY MODES) 


MODE 


WNHRMWDNHMPWDNHHMPWNHMWD MPWDNPWNPMWDHMWD 


REAL 


0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 


TABLE 1 (CONT) 


COEF. 


IMAG. 


0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
0.014860 
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SYMBOLS AND UNITS 


EX AXTAL COORDINATE OF FIELD POINT, 
FRACTION OF RADIUS 
ER RADIAL COORDINATE OF FIELD POINT, 


FRACTION OF RADIUS 
TETA ANGULAR COORDINATE OF FIELD POINT, 


RADIANS 
NT NONDIMENSIONAL TIME, NT=TIME*RPM/60 
KP PRESSURE COEFF., P/ (DENSITY* ( (RPM/ 


60) **2)* (DIAMETER **2) ) 


EX= -0.8000  ER= 1.2000 TETA= 0.0 
1 ~- TOTAL PRESSURE 
NT KP 
REAL IMAG 
0.0 0.109260 0.064822 
0.0625 ~0.002797 0.087227 
0.1250 ~0,101459 0.032849 
0.1875 ~0.086291 ~0.053374 
0.2500 0.115199 ~0.156751 
0.3125 0.085257 -~0.032765 
0.3750 0.021487 0.064308 
0.4375 ~0.051879 0.049902 
0.5000 ~0.0304070 ~-0.024493 
0.5625 ~0.045830 ~0.065644 
0.6250 0.036785 0.004191 
0.6875 0.028781 0.083445 
0.7500 -0,047159 0.067708 
0.8125 ~0.08048 ~0.061579 
0.8750 0.040268 ~0.055254 
0.9375 0.117604 0.020541 


TABLE 1 (cont) 
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2 - HARMONIC COMPONENTS 


HARM.NO. 


Oo ON WDB F WD NY FY OO 


~ Pe eH HY 
mF WwW DY F Oo 


REAL 
0.000563 
~0.000115 
0.008606 
0.069201 
-0.001695 
-0,011514 
~0.010340 
0.000716 
0.010076 
-0.001120 
-0,000015 
0.000433 
0.019964 
0.000831 
-0.000070 
0.018597 


KP 


IMAG 
0.000171 
0.001645 
0.000716 
0.047710 

-0,002692 
-0.011500 
0.001098 
0.001721 
0.001721 
-0.,009159 
0.000259 
0.000063 
-0.009159 
-0,000075 
~0,000130 
0.007825 


TABLE 1 (cont.) 
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the wake after about 7 revoiutions was found negligible, which agrees 
with calculation by Kerwin [18] for the zeroth harmonic. For field 
points down stream the wake integration is carried down 7 revolutions 
beyond the axial position of the field point. 

It is shown in Table 2 how the number of chordwise panels and 
the angular grid spacing affect the value of the harmonic components 
of the pressure. 

Figures 9 and 10 show the pressure variation and harmonic 
components at the field point due to the zeroth and third harmonic 
of loading respectively. 

Since we used an arbitrary loading’ distribution we cannot assess 
the accuracy of the results. This will have to be done at a later 


date. 





No. of panels 


Grid spacing 


HARM .NO. 
REAL 


0.0 


so WW W OO 


No. of panels 


Grid spacing 


HARM .NO. 
REAL 
0 0.0 
3 
6 -0.010340 
a 0.001120 


0.081273 


Boe 


A 


15 


4k degrees 


K 
P IMAG 


0.0 
0.054870 
0.000345 


0.002288 -0.017413 


20 


2 degrees 


K 
P IMAG 


0.0 


0.069201 0.047710 


0.001098 
-0.009159 


15 


2 degrees 


Kp 
REAL IMAG 


O70 0,0 


0.076421 0.053042 
-0.009018 0.000736 


0.001631 -0.012273 


TABLE 2 - EFFECT OF NUMBER OF CHORDWISE PANELS 


AND ANGULAR GRID SPACING 
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Field point: x/R 


r/R 


HARM.NO. 


© 


o DN WBN fF W DY & 


bec. Job —- PeeooURe AP UPSTREAM FIELD POINT DUE TO ZEROTH 


HARMONIC OF LOADING - HARMONIC COMPONENTS 


0.8 
1.2 


O. 


HARMONIC COMPONENTS 


real 

-0.059293 
-0.000003 
-0.000002 
~0.001161 
-~0.000003 

0.000000 
~0.000897 
-~0.000006 

0.000000 
-~0.000006 


imag. 

0.0 
-~0.000004 
~0.,000002 
~0.064603 
~0.000001 
~0.000001 
-0.000608 

0.000004 

0.000000 
-0.000004 
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with. 


Field point: x/R = 0.8 
BR = 1.2 
O= 0. 
-HARMONIC COMPONENTS 
HARM .NO, KP 
real imag. 
0 0.000563 0.000171 
1 -0.000115 0.001645 
2 0.008606 0.000716 
3 0.069201 0.047710 
dy -0.001695 ~0.002692 
5 ~0.001514 -0.001509 
6 ~0.010340 0.001098 
7 0.000716 0.001721 
8 0.000076 -0.000898 
-~0.001120 -0.009159 


FIG.10.b - PRESSURE AT UPSTREAM FIELD POINT DUE TO 
THIRD HARMONIC OF LOADING - HARMONIC 


COMPONENTS 
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Chapter V 
CONCLUSIONS AND RECOMMENDATIONS 


A computerized procedure has been developed for the calculation 
of the pressure field generated by a propeller operating in an un- 
steady flow. We Sigieve that it represents an important improvement 
over the existing methods based on lifting line theory and circum- 
ferentially uniform wake. It is also a useful complement to the 
recently published work by Jacobs et al. [19] which develops a 
different approach to the same problem. 

A considerable amount of time was put into the optimization of 
the program but it is probable that it can still be ameliorated. 

One of the most time consuming operations -- the integration of the 
streamwise vorticity in the wake -- was successfully reduced by the 
use of the Euler Transformation. The two other operations whose 
improvement would have an important effect on the over-all i mpeouenant 
of the program are the integration of the radial vorticity in the 

wake and the calculation of the velocity potential. 

Sample calculations were made to determine the sensibility of 
the results to the number of chordwise panels, angular grid spacing 
and how far downstream integrations were carried. Assessment of the 
accuracy of the results was not possible since we used an arbitrary 
loading distribution. 

Certainly, the next thing that needs to be done is to compare 


results with experimental data or with the results obtained in [19]. 





— a 


This requires that the loading distribution be obtained by any 
of the existing methods referred to before. The program accepts 


the input loading distribution in terms of pressure or circulation. 
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Appendix A 
VELOCITY POTENTIAL INDUCED BY A SINGLE VORTEX 


In this appendix we compare the velocity potential induced by 
a single reentrant vortex with the one induced by a concentrated 
doublet of equivalent strength. 


The velocity potential due to a vortex is given by 





S ih cos ¥ 
O = + = dA (A-1) 
a1 J, 52 


where 
S = vortex strength 
A = area enclosed by the vortex 
Y = angle between the normal to A and the line from A 
to the field point 
D = distance from A ee the field point. 

Equation (A-1) is similar to the one that gives the velocity 
potential induced by a uniform distribution over A of source doublets 
of strength S. It is clear that the strength of the equivalent con- 
centrated doublet must be equal to S x A. 

In order to simplify the calculations, we will consider a 
rectangular vortex centered at the origin, as shown in Fig. A-l. 


The velocity potential at p is found using eq. (A-1) 
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FIG. A-1 - ORIENTATION OF REENTRANT VORTEX 


The above integral can be evaluated analytically by integrating 
dike tly in one of the variables -- 6, say -- and then making the 


following changes of variables before performing the next integration: 


. / nea 

= a 

Reeeaye fe |" 

ps 2a 

u=vet(v +1)? 
2 
t = u 


The final result is 


_ 3 ln 1 p lias 1 : 
oe = aq bE Cean A - tan’ B) (tan’C - tan’ D) | (A-2) 
where 
De 
he Gia) 42° - (a) eee eee 
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ae sy ee ya Cie OG ee 
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z(x-a) 
D = (yta)"+2° - (yta)[(x-a) 4 (yta) 42°]? 
| 7 z(x - a) 


The velocity potential induced at p by a concentrated doublet 


2 
of strength 4a S, located at the origin is given by 


cae 
- a a ae, (A-3) 
W(x +y +2) 


tuto 


It is physically reasonable to predict that the largest dif- 
ferences between the values given by eqs. (A-1) and (A-2) will be found 
at the field points lying on the z axis. Figure A-2 shows how those 
differences vary with the distance along the z axis. 

It is apparent from the analysis of Fig. A-2 that the ratio of 
the two values approaches very rapidly the unity. For 2/a = 8 (i.e., 
for a distance 4 times greater than the side of the original vortex) 


the difference between the two values is less than one percent. 
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Appendix B 
COMPUTER PROGRAM 


1. General 

The program presented in this appendix is essentially the 
translation into FORTRAN IV of the numerical solution outlined in 
Chapter III. 

The initial section of the program in which the blade and wake 
are divided into panels is taken from Professor Kerwin's program, [18]. 
Also, the scheme used for the zeroth harmonic of the loading is based 
on Professor Kerwin's work. 

The program is written for use on the IBM System/360 digital 
computer. For different processors slight changes might have to be 
made. The amount of memory space required is about 150K bytes and 
the computation time for a single field point, including harmonic 
analysis of the pressure is about 4 minutes for the zeroth Seaseie 
and § minutes for higher harmonics, using a FORTRAN G Compiler. 

The complex harmonic analysis is performed by the subroutine 
HARM described in [23]. The final output variables are subscripted. 


This was done in anticipation of a plotting subroutine which can be 


added in the future. 
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2. List of Key Variables 


The following is a list of the principal variables in the program, 


arranged 


A(M,N) 


BOUNDA 


BOUNDT 


BWAKEA 


BWAKET 


CHOR (M) 


DHUB 


DIAM 


ER 


EX 


G(M) 


GB(M,N) 


GL (M,N) 


GLR(M) 


in alphabetical order. 


Coefficients of the Birnbaum modes of the loading function 
Axial velocity of the field point, induced by the bound 
vortices 

Tangential velocity at the field point, induced by the 
bound vortices 

Axial velocity at the field point, induced by the radial 
vortices in the wake 

Tangentical velocity at the field point, induced by the 
radial vortices in the wake 

Local chord at radius R(M) 

Hub diameter 

Propeller diameter 

Radial position of the field point 

Axial position of the field point 

Chordwise integral of the bound circulation at radius R(M) 
Bound circulation assumed concentrated on the N'th grid 
line in the M'th chordwise panel 

Wake circulation assumed concentrated between the 

KWAKE(M) + N and KWAKE(M) +N + 1 grid lines in the M'th chord- 
wise panel, for computing PTWK 


Advance coefficient at radius R(M) 





GLS 


Pp 


GR 


GT (M,N) 


IDENT 


IGB 


KD(N) 


KN (M) 


KWAKE (M) 


MN 


NBLADE 


NHARM 


NL 


NLE (M) 


NR 


NT 


NTE (CM) 


Seb 


Advance coefficient based on VS 


Circulation in the blade, assumed concentrated in the 
middle of each element of the lattice for computing PTBL 
Radial circulation in the wake 
Streamwise circulation in the blade, between the N and 
N+ 1 grid lines in the M'th panel 
Identification of the propeller 

Parameter that describes the loading function. If IGB=0 
the loading is in terms of bound circulation, if IGB > 0 
Bessie terms » the pressure jump in the blade 

Angle between the N'th blade and the YY axis 

Reduced frequency at radius R(M) 

Identification of the first grid line of the wake ar 
radius RZ(M) 

Number of equal spaces in the division of the span into 
chordwise panels 

Number of blades 

Harmonic number 

Largest integer multiple of the grid spacing with value 
less or equal to TMIN 

Identification of the first grid line in the M'th panel 
Number of loading modes 

Largest integer multiple of the grid spacing with value 
less or equal to TMAX 

Total number of grid lines in the blade 


°o 
Identification of the last grid in the blade in the M'th panel 





NTH 


NVV (M) 


PHN (M) 


PRS (M) 


Jase bE 


PTWK 


PTSC 


R(M) 


RPM 


RZ(M) 


SB(M,N) 


SL(M) 


SMA (N) 


SOURCA 


SOURCT 


tI 
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Parameter that defines the spacing between grid lines as 
NTH xX 2 degrees 

Total number of grid lines spanned by the blade in the 
M'th panel 

Phase angle of the NHARM'th harmonic of the inflow velocity 
Field point pressure when the first blade is at an angie 
of M*17//(8 NBLADE) with fixed vertical axis 

Time derivative of the velocity potential at the field 
point, induced by the circulation in the blade 

Time derivative of the velocity potential at the field 
point, induced by the circulation in the wake 

Time derivative of the velocity potential at the field 
point, induced by the source distribution 

Midradius of the M'th chordwise panel 

Rotation per minute 

Radius at the lower bound of the M'th panel 

Source strength concentrated on the N'th grid line in 
the M'th chordwise panel 

Leading edge at radius R(M) 

Derferential of the thickmess atethe N'th grid line, as 
percent of the max thickness in each panel 

Axial velocity at the field point, induced by the source 
distribution 

Tangential velocity at the field point, induced by the 


source distribution 





ST (M) 


STAR(M) 


eA 
TK(N) 


TL(M) 


TLE (M) 


TMIN 


TRAILA 


TRAILT 


TTE(M) 


TWAKEA 


TWAKET 


VFP 

VS 

VXO (M) 
XGL(N) 


XKN(N) 
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Trailing edge at radius R(M) 

Factor which multiplied by SMA(N) gives the source strength 
at the N'th grid line in the M'th chordwise panel 
Angular coordinate of the field point 

Standard thickness, fraction of the max thickness 
Amplitude of the streamwise vorticity in the wake at 
radius RZ(M) 

Angular coordinate of the leading edge at radius R(M) 
Extreme value of TTE(M) over all the panels 

Extreme value of TLE(M) over all the panels 

Axial velocity at the field point, induced by the 
streamwise vorticity in the blade 

Tangential velocity at the field Peat. induced by the 
en vorticity in the blade 

Angular coordinate of the trailing edge at radius R(M) 
Axial velocity at the field point, induced by the stream- 
wise vorticity in the wake 

Tangential velocity at the field point, induced by the 
streamwise vorticity in the wake 

Inflow velocity at the field point 

Ship speed 

Zeroth harmonic of the inflow velocity at radius R(M) 
Advance coefficient at radius XR(N) 


Reduced frequency at radius XR(N) 





X PHN (N) 


XR(N) 
XSL(N) 
XST(N) 


XSTAR(N) 


XTBETA(N) 
XTZ(N) 


XVXN (N) 
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Phase angle of the NHARM'th harmonic of the inflow 
velocity at radius XR(N) 

Input radius (see Input Description) 

Leading edge at radius XR(N) 

Trailing edge at radius XR(N) 

Factor of proportionality of the source strength at 
EPadwus XRCUN) 

Tangent of the advance angle at radius XR(N) 
Maximum thickness at radius XR(N) 


NHARM'th harmonic of the inflow velocity at radius XR(N) 
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co inpuc Desenmmptron 


The following sequence of input cards is required for the 


operation of the program: 


A IDENTIFICATION 1 card 

B GENERAL DESCRIPTION 1 card 

C BLADE OUTLINE 3 cards 

D WAKE CHARACTERISTICS 4. cards 

E LOADING DISTRIBUTION 21 cards max. 


F FIELD POINT DESCRIPTION 1 card 
A complete description of the above input cards is presented in 
Table B-l. The input radii referred to-in that table are schematically 


shown in the figure below. 


XR(1) XR(2) xR(9) XR(10) 


™..\ ee 


eee hs ee 


rp,/R 1.0 
DX 


FIG. Bel - INPUT RADII OF PROPELLER DATA 


DX is the Input Spacing, equal to % x interval from hub to tip. 
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